## ============================================================================
## 07_si_s2_audience_heterogeneity.R
##
## Purpose:   Treatment effect heterogeneity analysis for Study 1 using
##            Generalized Random Forests (causal forests). Estimates CATEs,
##            conducts omnibus tests for heterogeneity
##            Corresponds to SI Section S2.2.1.
##
## Inputs:    combined_data.rds (via 00_setup.R)
## Outputs:   grf_calibration_au_study1.rds,
##            grf_preds_au_study1.rds, 
##            grf_calibration_us_study1.rds,
##            grf_preds_us_study1.rds, 
##            Figures/fig_s13.pdf, 
##            Figures/fig_s14.pdf, 
##            Tables/table_s27.tex, 
##            Tables/table_s28.tex, 
## ============================================================================

source("00_setup.R")

#### ============== Estimate causal forests: Australia ============== ####
au_covs <-
  c(
    "x_gender_f",
    "x_age_cat",
    "x_state_comb",
    "x_edu_college",
    "x_native",
    "x_employ_cat",
    "x_hhi",
    "x_pid_comb2",
    "x_ideo_cat"
  )

au_outcomes <- c(
  "y_guilt_idx_s",
  "y_cbr_idx_s",
  "y_symbolic_idx_s",
  "y_substantive_idx_s",
  "y_voice_vote_s",
  "y_infoseek_s",
  "y_acknowledge_s"
)

# NB: GRF only accepts rows with complete cases on DVs and all
# au_covs.
est_dat <-
  combined_dat %>%
  filter(sample == "AU") %>%
  select(all_of(au_covs), Z_land_info, all_of(au_outcomes))

index <- complete.cases(est_dat)
X <- est_dat[index, c(au_covs)]

W <- as.numeric(est_dat$Z_land_info[index])

## Recode au_covs to numeric
dummies <- model.matrix(~ x_age_cat - 1, data = X)
X <- cbind(X, dummies)

dummies <- model.matrix(~ x_hhi - 1, data = X)
X <- cbind(X, dummies)

dummies <- model.matrix(~ x_ideo_cat - 1, data = X)
X <- cbind(X, dummies)

dummies <- model.matrix(~ x_pid_comb2 - 1, data = X)
X <- cbind(X, dummies)

dummies <- model.matrix(~ x_employ_cat - 1, data = X)
X <- cbind(X, dummies)

dummies <- model.matrix(~ x_state_comb - 1, data = X)
X <- cbind(X, dummies)

## Keep numeric au_covs for GRF:
nums <- unlist(lapply(X, is.numeric))
X <- as.matrix(X[, nums])

# Train a causal forest using the observed outcome (Y_Obs), the matrix
# of pre-treatment au_covs (X), and reatment assignment vector (W).
n <- nrow(est_dat[index,])

#### ---------- Support for the voice (AU only) ----
Y_Obs <- est_dat$y_voice_vote_s[index]
tau_voice <-
  causal_forest(
    Y = Y_Obs,
    X = X,
    W = W,
    W.hat = 0.5,
    honesty = TRUE,
    num.trees = 5000,
    seed = 123
  )

# Estimate the treatment effects for all units using out-of-bag prediction.
tau_hat_voice <- predict(tau_voice, estimate.variance = TRUE)

# Compute estimated SEs treatment effects for predictions.
sigma_hat_voice <- sqrt(tau_hat_voice$variance.estimates)

# Store the estimated treatment effects and CIs in a data frame.
grf_voice <-
  cbind(data.frame(
    tau_hat = tau_hat_voice$predictions,
    ci_upper = tau_hat_voice$predictions + 1.96 * sigma_hat_voice,
    ci_lower = tau_hat_voice$predictions - 1.96 * sigma_hat_voice),
    X
  ) %>%
  mutate(outcome = "Voting for Referendum")

#### ---------- Information seeking ----
Y_Obs <- est_dat$y_infoseek_s[index]
tau_infoseek <-
  causal_forest(
    Y = Y_Obs,
    X = X,
    W = W,
    W.hat = 0.5,
    honesty = TRUE,
    num.trees = 5000,
    seed = 123
  )

# Estimate the treatment effects for all units using out-of-bag prediction.
tau_hat_infoseek <- predict(tau_infoseek, estimate.variance = TRUE)

# Compute estimated SEs treatment effects for predictions.
sigma_hat_infoseek <- sqrt(tau_hat_infoseek$variance.estimates)

# Store the estimated treatment effects and CIs in a data frame.
grf_infoseek <-
  cbind(data.frame(
    tau_hat = tau_hat_infoseek$predictions,
    ci_upper = tau_hat_infoseek$predictions + 1.96 * sigma_hat_infoseek,
    ci_lower = tau_hat_infoseek$predictions - 1.96 * sigma_hat_infoseek),
    X
  ) %>%
  mutate(outcome = "Information Seeking")

#### ---------- Binary belief ----
Y_Obs <- est_dat$y_acknowledge_s[index]
tau_binbelief <-
  causal_forest(
    Y = Y_Obs,
    X = X,
    W = W,
    W.hat = 0.5,
    honesty = TRUE,
    num.trees = 5000,
    seed = 123
  )

# Estimate the treatment effects for all units using out-of-bag prediction.
tau_hat_binbelief <- predict(tau_binbelief, estimate.variance = TRUE)

# Compute estimated SEs treatment effects for predictions.
sigma_hat_binbelief <- sqrt(tau_hat_binbelief$variance.estimates)

# Store the estimated treatment effects and CIs in a data frame.
grf_binbelief <-
  cbind(data.frame(
    tau_hat = tau_hat_binbelief$predictions,
    ci_upper = tau_hat_binbelief$predictions + 1.96 * sigma_hat_binbelief,
    ci_lower = tau_hat_binbelief$predictions - 1.96 * sigma_hat_binbelief),
    X
  ) %>%
  mutate(outcome = "Indigenous Land")

#### ---------- Collective Guilt ----
Y_Obs <- est_dat$y_guilt_idx_s[index]
tau_guilt <-
  causal_forest(
    Y = Y_Obs,
    X = X,
    W = W,
    W.hat = 0.5,
    honesty = TRUE,
    num.trees = 5000,
    seed = 123
  )

# Estimate the treatment effects for all units using out-of-bag prediction.
tau_hat_guilt <- predict(tau_guilt, estimate.variance = TRUE)

# Compute estimated SEs treatment effects for predictions.
sigma_hat_guilt <- sqrt(tau_hat_guilt$variance.estimates)

# Store the estimated treatment effects and CIs in a data frame.
grf_guilt <-
  cbind(data.frame(
    tau_hat = tau_hat_guilt$predictions,
    ci_upper = tau_hat_guilt$predictions + 1.96 * sigma_hat_guilt,
    ci_lower = tau_hat_guilt$predictions - 1.96 * sigma_hat_guilt),
    X
  ) %>%
  mutate(outcome = "Collective Guilt")

#### ---------- Colorblind Racism ----
Y_Obs <- est_dat$y_cbr_idx_s[index]
tau_cbr <-
  causal_forest(
    Y = Y_Obs,
    X = X,
    W = W,
    W.hat = 0.5,
    honesty = TRUE,
    num.trees = 5000,
    seed = 123
  )

# Estimate the treatment effects for all units using out-of-bag prediction.
tau_hat_cbr <- predict(tau_cbr, estimate.variance = TRUE)

# Compute estimated SEs treatment effects for predictions.
sigma_hat_cbr <- sqrt(tau_hat_cbr$variance.estimates)

# Store the estimated treatment effects and CIs in a data frame.
grf_cbr <-
  cbind(data.frame(
    tau_hat = tau_hat_cbr$predictions,
    ci_upper = tau_hat_cbr$predictions + 1.96 * sigma_hat_cbr,
    ci_lower = tau_hat_cbr$predictions - 1.96 * sigma_hat_cbr),
    X
  ) %>%
  mutate(outcome = "Colorblind Racism")

#### ---------- Symbolic Reparations ----
Y_Obs <- est_dat$y_symbolic_idx_s[index]
tau_symrep <-
  causal_forest(
    Y = Y_Obs,
    X = X,
    W = W,
    W.hat = 0.5,
    honesty = TRUE,
    num.trees = 5000,
    seed = 123
  )

# Estimate the treatment effects for all units using out-of-bag prediction.
tau_hat_symrep <- predict(tau_symrep, estimate.variance = TRUE)

# Compute estimated SEs treatment effects for predictions.
sigma_hat_symrep <- sqrt(tau_hat_symrep$variance.estimates)

# Store the estimated treatment effects and CIs in a data frame.
grf_symrep <-
  cbind(data.frame(
    tau_hat = tau_hat_symrep$predictions,
    ci_upper = tau_hat_symrep$predictions + 1.96 * sigma_hat_symrep,
    ci_lower = tau_hat_symrep$predictions - 1.96 * sigma_hat_symrep),
    X
  ) %>%
  mutate(outcome = "Symbolic Reparations")

#### ---------- Substantive Reparations ----
Y_Obs <- est_dat$y_substantive_idx_s[index]
tau_subrep <-
  causal_forest(
    Y = Y_Obs,
    X = X,
    W = W,
    W.hat = 0.5,
    honesty = TRUE,
    num.trees = 5000,
    seed = 123
  )

# Estimate the treatment effects for all units using out-of-bag prediction.
tau_hat_subrep <- predict(tau_subrep, estimate.variance = TRUE)

# Compute estimated SEs treatment effects for predictions.
sigma_hat_subrep <- sqrt(tau_hat_subrep$variance.estimates)

# Store the estimated treatment effects and CIs in a data frame.
grf_subrep <-
  cbind(data.frame(
    tau_hat = tau_hat_subrep$predictions,
    ci_upper = tau_hat_subrep$predictions + 1.96 * sigma_hat_subrep,
    ci_lower = tau_hat_subrep$predictions - 1.96 * sigma_hat_subrep),
    X
  ) %>%
  mutate(outcome = "Substantive Reparations")

#### ------ Stack results and export ----

au_forest_test_df <-
  bind_rows(
    "Indigenous Land" = broom::tidy(test_calibration(tau_binbelief)),
    "Information Seeking" = broom::tidy(test_calibration(tau_infoseek)),
    "Voting for Referendum" = broom::tidy(test_calibration(tau_voice)),
    "Collective Guilt" = broom::tidy(test_calibration(tau_guilt)),
    "Colorblind Racism" = broom::tidy(test_calibration(tau_cbr)),
    "Symbolic Reparations" = broom::tidy(test_calibration(tau_symrep)),
    "Substantive Reparations" = broom::tidy(test_calibration(tau_subrep)),
    .id = "outcome"
  ) %>%
  mutate(outcome = factor(
    outcome,
    levels = c(
      "Voting for Referendum",
      "Substantive Reparations",
      "Symbolic Reparations",
      "Colorblind Racism",
      "Collective Guilt",
      "Information Seeking",
      "Indigenous Land"
    ),
    labels = c(
      "Voice Referendum*",
      "Substantive Reparations",
      "Symbolic Reparations",
      "Colorblind Racism",
      "Collective Guilt",
      "Information Seeking",
      "Indigenous Land"
    )
  )) %>%
  mutate(across(where(is.numeric), \(x) round(x, 2)))

au_forest_test_df %>% 
  filter(term == "differential.forest.prediction") %>% 
  arrange(desc(outcome))

saveRDS(au_forest_test_df, "grf_calibration_au_study1.rds")

# Order treatment effects by size
grf_preds_au <-
  bind_rows(
    grf_binbelief,
    grf_voice,
    grf_infoseek,
    grf_guilt,
    grf_cbr,
    grf_symrep,
    grf_subrep
  ) %>%
  mutate(outcome = factor(
    outcome,
    levels = c(
      "Voting for Referendum",
      "Substantive Reparations",
      "Symbolic Reparations",
      "Colorblind Racism",
      "Collective Guilt",
      "Information Seeking",
      "Indigenous Land"
    ),
    labels = c(
      "Voice Referendum*",
      "Substantive Reparations",
      "Symbolic Reparations",
      "Colorblind Racism",
      "Collective Guilt",
      "Information Seeking",
      "Indigenous Land"
    )
  )) %>%
  group_by(outcome) %>%
  mutate(order = rank(tau_hat, ties.method = "first")) %>%
  ungroup()

saveRDS(grf_preds_au, "grf_preds_au_study1.rds")

#### ============== Generate Fig. S13, Table S27 ============== ####

grf_preds_au <- read_rds("grf_preds_au_study1.rds")
grf_calibrate_au <- read_rds("grf_calibration_au_study1.rds")

## Plot treatment effects 
pred_au <-
  grf_preds_au %>%
  mutate(
    outcome = factor(outcome,
                     levels = c("Indigenous Land", "Information Seeking",
                                "Collective Guilt", "Colorblind Racism",
                                "Symbolic Reparations",
                                "Substantive Reparations",
                                "Voice Referendum*"),
                     labels = c("Indigenous Land", "Information Seeking",
                                "Collective Guilt", "Colorblind Racism",
                                "Symbolic Reparations",
                                "Substantive Reparations",
                                "Voice Referendum"))
  ) %>%
  ggplot(., aes(x = tau_hat, y = order)) +
  facet_wrap( ~ outcome) +
  geom_errorbar(aes(xmin = ci_lower, xmax = ci_upper),
                col = "dimgrey",
                alpha = 0.2,
                width = 0,
                linewidth = 0.5,
                orientation = "y") +
  geom_point(size = .5, pch = 20,
             alpha = 0.8,
             fill = "black",
             col = "black"
  ) +
  geom_vline(xintercept = 0, linewidth = 0.5, lty = 2) +
  labs(title = "",
       x = "",
       y = "") +
  scale_y_continuous(expand = c(0.01,0.01)) +
  scale_x_continuous(expand = c(0.05,0.05)) +
  theme_tufte() +
  theme(axis.ticks.y = element_blank(),
        axis.ticks.x = element_line(linewidth = 0.25),
        strip.text = element_text(size = 11),
        axis.line.x = element_line(linewidth = 0.25),
        axis.title.x = element_text(size = 10),
        axis.text.y = element_blank(),
        axis.line.y = element_blank(),
        panel.border = element_rect(fill = NA),
        strip.background = element_rect(fill = "lightgrey"),
        plot.margin = unit(c(t = -0.5, r = 0.5, b = -0.5, l = -0.5), "lines")
  )

ggsave(
  pred_au,
  filename = paste0(fig_dir, "fig_s13.pdf"),
  width = 7,
  height = 10
)

## Export table of results from omnibus tests
grf_calibrate_au %>%
  filter(term == "differential.forest.prediction") %>%
  mutate(p.adj = p.adjust(p.value, method = "bonferroni")) %>%
  arrange(desc(outcome)) %>%
  select(-term) %>%
  xtable(.,) %>%
  print(
    include.rownames = FALSE,
    include.colnames = FALSE,
    hline.after = c(),
    only.contents = TRUE,
    type = "latex",
    file = paste0(tab_dir, "table_s27.tex")
  )

#### ============== Estimate causal forests: United States ============== ####
us_covs <-
  c(
    "x_gender_f",
    "x_age_cat",
    "x_ethnicity",
    "x_region",
    "x_edu_college",
    "x_native",
    "x_employ_cat",
    "x_hhi",
    "x_pid_comb2"
  )

us_outcomes <- c(
  "y_guilt_idx_s",
  "y_cbr_idx_s",
  "y_symbolic_idx_s",
  "y_substantive_idx_s",
  "y_infoseek_s",
  "y_acknowledge_s"
)

# NB: GRF only accepts rows with complete cases on DVs and all
# us_covs.
est_dat <-
  combined_dat %>%
  filter(sample == "US") %>%
  select(all_of(us_covs), Z_land_info, all_of(us_outcomes))

index <- complete.cases(est_dat)
X <- est_dat[index, c(us_covs)]

W <- as.numeric(est_dat$Z_land_info[index])

## Recode us_covs to numeric
dummies <- model.matrix(~ x_age_cat - 1, data = X)
X <- cbind(X, dummies)

dummies <- model.matrix(~ x_hhi - 1, data = X)
X <- cbind(X, dummies)

dummies <- model.matrix(~ x_ethnicity - 1, data = X)
X <- cbind(X, dummies)

dummies <- model.matrix(~ x_region - 1, data = X)
X <- cbind(X, dummies)

dummies <- model.matrix(~ x_pid_comb2 - 1, data = X)
X <- cbind(X, dummies)

dummies <- model.matrix(~ x_employ_cat - 1, data = X)
X <- cbind(X, dummies)

## Keep numeric us_covs for GRF:
nums <- unlist(lapply(X, is.numeric))
X <- as.matrix(X[, nums])

# Step 1. Train a causal forest using the observed outcome (Y_Obs), the matrix
# of pre-treatment us_covs (X), and reatment assignment vector (W).
set.seed(123)
n <- nrow(est_dat[index,])

#### ---------- Information seeking ----
Y_Obs <- est_dat$y_infoseek_s[index]
tau_infoseek <-
  causal_forest(
    Y = Y_Obs,
    X = X,
    W = W,
    W.hat = 0.5,
    honesty = TRUE,
    num.trees = 5000,
    seed = 123
  )

# Estimate the treatment effects for all units using out-of-bag prediction.
tau_hat_infoseek <- predict(tau_infoseek, estimate.variance = TRUE)

# Compute estimated SEs treatment effects for predictions.
sigma_hat_infoseek <- sqrt(tau_hat_infoseek$variance.estimates)

# Store the estimated treatment effects and CIs in a data frame.
grf_infoseek <-
  cbind(data.frame(
    tau_hat = tau_hat_infoseek$predictions,
    ci_upper = tau_hat_infoseek$predictions + 1.96 * sigma_hat_infoseek,
    ci_lower = tau_hat_infoseek$predictions - 1.96 * sigma_hat_infoseek),
    X
  ) %>%
  mutate(outcome = "Information Seeking")

#### ---------- Binary belief ----
Y_Obs <- est_dat$y_acknowledge_s[index]
tau_binbelief <-
  causal_forest(
    Y = Y_Obs,
    X = X,
    W = W,
    W.hat = 0.5,
    honesty = TRUE,
    num.trees = 5000,
    seed = 123
  )

# Estimate the treatment effects for all units using out-of-bag prediction.
tau_hat_binbelief <- predict(tau_binbelief, estimate.variance = TRUE)

# Compute estimated SEs treatment effects for predictions.
sigma_hat_binbelief <- sqrt(tau_hat_binbelief$variance.estimates)

# Store the estimated treatment effects and CIs in a data frame.
grf_binbelief <-
  cbind(data.frame(
    tau_hat = tau_hat_binbelief$predictions,
    ci_upper = tau_hat_binbelief$predictions + 1.96 * sigma_hat_binbelief,
    ci_lower = tau_hat_binbelief$predictions - 1.96 * sigma_hat_binbelief),
    X
  ) %>%
  mutate(outcome = "Indigenous Land")

#### ---------- Collective Guilt ----
Y_Obs <- est_dat$y_guilt_idx_s[index]
tau_guilt <-
  causal_forest(
    Y = Y_Obs,
    X = X,
    W = W,
    W.hat = 0.5,
    honesty = TRUE,
    num.trees = 5000,
    seed = 123
  )

# Estimate the treatment effects for all units using out-of-bag prediction.
tau_hat_guilt <- predict(tau_guilt, estimate.variance = TRUE)

# Compute estimated SEs treatment effects for predictions.
sigma_hat_guilt <- sqrt(tau_hat_guilt$variance.estimates)

# Store the estimated treatment effects and CIs in a data frame.
grf_guilt <-
  cbind(data.frame(
    tau_hat = tau_hat_guilt$predictions,
    ci_upper = tau_hat_guilt$predictions + 1.96 * sigma_hat_guilt,
    ci_lower = tau_hat_guilt$predictions - 1.96 * sigma_hat_guilt),
    X
  ) %>%
  mutate(outcome = "Collective Guilt")

#### ---------- Colorblind Racism ----
Y_Obs <- est_dat$y_cbr_idx_s[index]
tau_cbr <-
  causal_forest(
    Y = Y_Obs,
    X = X,
    W = W,
    W.hat = 0.5,
    honesty = TRUE,
    num.trees = 5000,
    seed = 123
  )

# Estimate the treatment effects for all units using out-of-bag prediction.
tau_hat_cbr <- predict(tau_cbr, estimate.variance = TRUE)

# Compute estimated SEs treatment effects for predictions.
sigma_hat_cbr <- sqrt(tau_hat_cbr$variance.estimates)

# Store the estimated treatment effects and CIs in a data frame.
grf_cbr <-
  cbind(data.frame(
    tau_hat = tau_hat_cbr$predictions,
    ci_upper = tau_hat_cbr$predictions + 1.96 * sigma_hat_cbr,
    ci_lower = tau_hat_cbr$predictions - 1.96 * sigma_hat_cbr),
    X
  ) %>%
  mutate(outcome = "Colorblind Racism")

#### ---------- Symbolic Reparations ----
Y_Obs <- est_dat$y_symbolic_idx_s[index]
tau_symrep <-
  causal_forest(
    Y = Y_Obs,
    X = X,
    W = W,
    W.hat = 0.5,
    honesty = TRUE,
    num.trees = 5000,
    seed = 123
  )

# Estimate the treatment effects for all units using out-of-bag prediction.
tau_hat_symrep <- predict(tau_symrep, estimate.variance = TRUE)

# Compute estimated SEs treatment effects for predictions.
sigma_hat_symrep <- sqrt(tau_hat_symrep$variance.estimates)

# Store the estimated treatment effects and CIs in a data frame.
grf_symrep <-
  cbind(data.frame(
    tau_hat = tau_hat_symrep$predictions,
    ci_upper = tau_hat_symrep$predictions + 1.96 * sigma_hat_symrep,
    ci_lower = tau_hat_symrep$predictions - 1.96 * sigma_hat_symrep),
    X
  ) %>%
  mutate(outcome = "Symbolic Reparations")
#### ---------- Substantive Reparations ----
Y_Obs <- est_dat$y_substantive_idx_s[index]
tau_subrep <-
  causal_forest(
    Y = Y_Obs,
    X = X,
    W = W,
    W.hat = 0.5,
    honesty = TRUE,
    num.trees = 5000,
    seed = 123
  )

# Estimate the treatment effects for all units using out-of-bag prediction.
tau_hat_subrep <- predict(tau_subrep, estimate.variance = TRUE)

# Compute estimated SEs treatment effects for predictions.
sigma_hat_subrep <- sqrt(tau_hat_subrep$variance.estimates)

# Store the estimated treatment effects and CIs in a data frame.
grf_subrep <-
  cbind(data.frame(
    tau_hat = tau_hat_subrep$predictions,
    ci_upper = tau_hat_subrep$predictions + 1.96 * sigma_hat_subrep,
    ci_lower = tau_hat_subrep$predictions - 1.96 * sigma_hat_subrep),
    X
  ) %>%
  mutate(outcome = "Substantive Reparations")

#### ------ Stack results and export ----

us_forest_test_df <-
  bind_rows(
    "Indigenous Land" = broom::tidy(test_calibration(tau_binbelief)),
    "Information Seeking" = broom::tidy(test_calibration(tau_infoseek)),
    "Collective Guilt" = broom::tidy(test_calibration(tau_guilt)),
    "Colorblind Racism" = broom::tidy(test_calibration(tau_cbr)),
    "Symbolic Reparations" = broom::tidy(test_calibration(tau_symrep)),
    "Substantive Reparations" = broom::tidy(test_calibration(tau_subrep)),
    .id = "outcome"
  ) %>%
  mutate(outcome = factor(
    outcome,
    levels = c(
      "Indigenous Land",
      "Information Seeking",
      "Collective Guilt",
      "Colorblind Racism",
      "Symbolic Reparations",
      "Substantive Reparations"
    )
  )) %>%
  mutate(across(where(is.numeric), \(x) round(x, 2)))

us_forest_test_df %>% 
  filter(term == "differential.forest.prediction") %>% 
  arrange(outcome)

saveRDS(us_forest_test_df, "grf_calibration_us_study1.rds")

# Order treatment effects by size
grf_preds_us <-
  bind_rows(
    grf_binbelief,
    grf_infoseek,
    grf_guilt,
    grf_cbr,
    grf_symrep,
    grf_subrep
  ) %>%
  mutate(outcome = factor(
    outcome,
    levels = c(
      "Indigenous Land",
      "Information Seeking",
      "Collective Guilt",
      "Colorblind Racism",
      "Symbolic Reparations",
      "Substantive Reparations"
    )
  )) %>%
  group_by(outcome) %>%
  mutate(order = rank(tau_hat, ties.method = "first")) %>%
  ungroup()

saveRDS(grf_preds_us, "grf_preds_us_study1.rds")

#### ============== Generate Fig. S14, Table S28 ============== ####
grf_preds_us <- read_rds("grf_preds_us_study1.rds")
grf_calibrate_us <- read_rds("grf_calibration_us_study1.rds")

## Plot treatment effects 
pred_us <-
  grf_preds_us %>%
  ggplot(., aes(x = tau_hat, y = order)) +
  facet_wrap( ~ outcome) +
  geom_errorbar(aes(xmin = ci_lower, xmax = ci_upper),
                 col = "dimgrey",
                 alpha = 0.2,
                 width = 0,
                 linewidth = 0.5,
                 orientation = "y") +
  geom_point(size = .5, pch = 20,
             alpha = 0.8,
             fill = "black",
             col = "black"
  ) +
  geom_vline(xintercept = 0, linewidth = 0.5, lty = 2) +
  labs(title = "",
       x = "",
       y = "") +
  scale_y_continuous(expand = c(0.01,0.01)) +
  scale_x_continuous(expand = c(0.05,0.05)) +
  theme_tufte() +
  theme(axis.ticks.y = element_blank(),
        axis.ticks.x = element_line(linewidth = 0.25),
        strip.text = element_text(size = 11),
        axis.line.x = element_line(linewidth = 0.25),
        axis.title.x = element_text(size = 10),
        axis.text.y = element_blank(),
        axis.line.y = element_blank(),
        panel.border = element_rect(fill = NA),
        strip.background = element_rect(fill = "lightgrey"),
        plot.margin = unit(c(t = -0.5, r = 0.5, b = -0.5, l = -0.5), "lines")
  )

ggsave(
  pred_us,
  filename = paste0(fig_dir, "fig_s14.pdf"),
  width = 7,
  height = 7
)

## Table of results from omnibus tests
grf_calibrate_us %>%
  filter(!str_detect(term, "mean.forest")) %>%
  mutate(p.adj = p.adjust(p.value, method = "bonferroni")) %>%
  arrange(outcome) %>%
  select(-term) %>%
  xtable(.,) %>%
  print(
    include.rownames = FALSE,
    include.colnames = FALSE,
    hline.after = c(),
    only.contents = TRUE,
    type = "latex",
    file = paste0(tab_dir, "table_s28.tex")
  )


## ---- Session Info -----------------------------------------------------------
sessionInfo()